home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YU.C < prev    next >
Text File  |  1988-12-06  |  19KB  |  674 lines
  1. /******************* (C) 1986,7,8 by JAMES A. YORKE **************************/
  2. /**********************************  YU.C  ***********************************/
  3.  
  4.  
  5. #include "yinclud.h"
  6.  
  7. static double   diagonal;    /* for unstable manifold routines in this file 
  8.                 */
  9. /* int goodIterate; the number of iterates before leaving the box 
  10.         (set in terms of diameters[] ) */
  11. static double   ddx,
  12.                 ddy;
  13. static double   PreviousAdd,
  14.                 oneOverX,
  15.                 oneOverY;
  16. static int      multiplier;
  17. static double   jump,
  18.                 oldX,
  19.                 oldY,
  20.                 dif;
  21. static int      itermax;
  22.  
  23. LineIterator(iterMax, iterMultiplier)/* The images of the line from  
  24.                    x1_line,y1_line   to   x2_line,y2_line are
  25.                    drawn, where iterMultiplier is a number >=1
  26.                    such that the line is iterated in multiples
  27.                    of iterMultiplier; in the computation of
  28.                    unstable manifolds of periodic orbits, it
  29.                    should be 2*period where period is the
  30.                    period of the orbit; the factor of two is
  31.                    because the orbit might have a negative
  32.                    expanding eigenvalue; iterMax is the maximum
  33.                    number of times "iter" is increased; the
  34.                    actual maximum number of iterates is
  35.                    iterMax*iterMultiplier; routines in this
  36.                    file: set_line(),LineForB(), and LineForU()
  37.                    set the points(line_xi,line_yi) for i = 1
  38.                    and 2.             */
  39. int     iterMax,
  40.         iterMultiplier;
  41. {
  42.     int     isInBox(), level0;
  43.     int     checkFreq();
  44.     double  minAdd = 1;
  45.     double  getAdd();
  46.     double  fracY[EQUATIONS];
  47.     double  oldLength = 0.,
  48.             entropy,
  49.             difX,
  50.             difY;
  51.  
  52.  /* the following are for restarting this routine */
  53.  /* dot is just for keeping track of how much has been done */
  54.  
  55.     if(refreshFlag == YES) {
  56.         if(dot > 0 && dot < dots)
  57.                 /* dot >0 is most likely to occur when using
  58.                    command "cont" */
  59.             ;
  60.  
  61.     }
  62.     else {
  63.         dot = 0;
  64.         iter = 1;
  65.     }
  66.  
  67.  /* dot is just for keeping track of how much has been done */
  68.     initTime();
  69.     itermax = iterMax;
  70.     multiplier = iterMultiplier;/* we need this number globally */
  71.     level0 = ++level;
  72.     if(level == PROCESS && refreshFlag == NO)
  73.         boot_crt();
  74.     store(y, y + eqn0, eqn1);
  75.     oldX = y[X_coord];
  76.     oldY = y[Y_coord];
  77.     jump = USTEP;        /* this is a fraction of the screen's diagonal
  78.                    length(currently .001) indicating how big a
  79.                    jump is allowed so that we do not move more
  80.                    than 1 pixel */
  81.  
  82.     difX = line_x1 - line_x2;
  83.     difY = line_y1 - line_y2;
  84.     length = sqrt(difX * difX + difY * difY);
  85.                 /* this will tally the total length plotted so
  86.                    far; initially it is the length of the
  87.                    initial line segment */
  88.     diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
  89.             + (Y_upper - Y_lower) * (Y_upper - Y_lower));
  90.     boxConstants();    /* this defines constants needed by isInBox()
  91.                    for checking if y[] is in the box */
  92.     oldGoodIt = 0;
  93.     goodIterate = 1;
  94.     iter = 1;
  95.     add = USTEP;
  96.     store(fracY, y, dim);    /* needed for getAdd() */
  97.     add = getAdd(add, fracY);/* one iterate; only place called */
  98.     for(; level == level0 && iter <= itermax; iter++) {
  99.         goodIterate = iter;
  100.     /* first we compute the initial jump size "add" */
  101.         if(checkFreq () == YES) {
  102.                 /* means terminate; checkFreq also handles
  103.                    occasional transmissions of picture to disk 
  104.                 */
  105.             PRINT "process termination: OUT OF TIME \n");
  106.             level = level0 - 1;
  107.             continue;
  108.         }
  109.         scr_rowcol(4, 0);
  110.         if(printer >= 2) {
  111.             if(iter > 1 && length > 0.&& oldLength > 0.) {
  112.                 entropy = log(length / oldLength);
  113.                 PRINT
  114.                     "length(%d)=length of %dth image of segment=%lf   \n",
  115.                     iter - 1, iter - 1, length * diagonal);
  116.                 PRINT
  117.                     "Crude Estimate of topological entropy = log(length(%d)/length(%d))=%5.5lf  \n"
  118.                     ,iter - 1, iter - 2, entropy);
  119.             }
  120.             else
  121.                 PRINT "length(0)=%lf \n", length);
  122.  
  123.         }
  124.         oldLength = length;
  125.         length = 0.;
  126.  
  127.     /* MAIN INNER LOOP */
  128.         if(refreshFlag == YES)
  129.             refreshFlag = NO;
  130.         else
  131.             frac = 0;/* this occurs always except possibly once */
  132.         for(; level == level0 && frac < 1.0;
  133.                 frac = ((frac + add >= 1.0) ? 1.0 : frac + add)) {
  134.             checkForInterrupt();
  135.             dot++;    /* can be deleted; just for keeping track of
  136.                    how much has been done and how fast it is */
  137.         /* when frac == 0, goodIterate and oldGoodIterate will not be
  138.            equal since goodIterate has just been increased beyond
  139.            previous values */
  140.  
  141.             getGoodIt();
  142.                 /* this is the only place in the "for(frac..."
  143.                    loop that goodIterate is set; this also sets
  144.                    y according to frac and iterates as many
  145.                    times up to iter as possible while staying
  146.                    in the box; at the end, y[] is this final
  147.                    image point that is in the box and
  148.                    goodIterate is the number of good iterates 
  149.                 */
  150.  
  151.             while(goodIterate >= oldGoodIt + 2) {
  152.                 scr_rowcol(7, 0);
  153.                 PRINT
  154.                     "good Iterate just jumped by 2\n");
  155.                 frac = frac - add;
  156.                 add = add / 5.;
  157.                 frac = frac + add;
  158.                 getGoodIt();
  159.             }
  160.             store(fracY, y, dim);/* fracY is used by revise... */
  161.  
  162.  
  163.             if (printer >= 3)
  164.             {
  165.                 scr_rowcol(9,0);
  166.                 PRINT 
  167. "%d good; iter=%d frac=%lf,add=%le  \n",goodIterate,iter,frac,add);
  168.             }
  169.             if (printer >= 1 && minAdd > add)
  170.             {
  171.                 scr_rowcol(8,0);
  172.                 PRINT 
  173. "min add so far =%le             \n",add);
  174.             }
  175.             minAdd = (minAdd > add ? add : minAdd ) ;
  176.  
  177.             reviseAddAndJump(goodIterate, fracY);
  178.  
  179.         /* if(goodIterate <= iter-1) we test to see if frac+add gives
  180.            a point that can be iterated one more time; if it can be,
  181.            then we will be decrease add by implementing the instruction
  182.            : reviseAddAndJump(goodIterate+1,fracY) which should produce
  183.            a smaller "add" and just to be safer, we "add = add/3.;" */
  184.  
  185.             if(goodIterate <= iter - 1) {
  186.                 process(1);
  187.                 /* one more iterate from y[] as finally set by
  188.                    revise...() */
  189.                 if(isInBox (y[X_coord], y[Y_coord]) == YES) {
  190.                     store(y, fracY, dim);
  191.                     process(1);
  192.                 /* one more iterate from y[] as finally set by
  193.                    GetGoodIt() */
  194.                     store(fracY, y, dim);
  195.                 /* now fracY is as set by GetGoodIt() PLUS 1
  196.                    iterate */
  197.                     reviseAddAndJump(goodIterate + 1, fracY);
  198.                     add = add / 3.;
  199.                 }
  200.             }
  201.             if(goodIterate == iter || oldGoodIt == iter) {
  202.                 store(y, fracY, dim);
  203.                 plotSeg();
  204.                 length = length + dif;
  205.             }
  206.  
  207.             oldGoodIt = goodIterate;
  208.         /* the following two coordinates result from iterating y[] iter
  209.            times PROVIDED goodIterate = either iter or iter -1; then
  210.            they can be used for plotting */
  211.             oldX = y[X_coord];
  212.             oldY = y[Y_coord];
  213.         }
  214.     }
  215.     level = level0 - 1;
  216.     if(level < PROCESS)
  217.         MainMenu();
  218.  
  219. }
  220.  
  221.  
  222. double  jumpSize() {        /* returns the maximum of the horizontal and
  223.                    vertical fractions of the screen between y
  224.                    and(oldX,oldY) */
  225.     double  fabs(), difX, difY;
  226.  
  227.     difX = fabs((y[X_coord] - oldX) * oneOverX);
  228.     difY = fabs((y[Y_coord] - oldY) * oneOverY);
  229.     if(difX > difY)
  230.         return(difX);
  231.     return(difY);
  232. }
  233.  
  234. reviseAddAndJump(its, pY)    /* given an estimate "add" of how much frac
  235.                    should be changed by, try it out using
  236.                    MakeAJump() which tells us what fraction of
  237.                    the screen we actually move; if that move is
  238.                    too big, (that is, if the updated estimate
  239.                    is less than .8 of add), then do it again
  240.                    using the new add; the desired jump is
  241.                    jump*numPixels so we change add by the
  242.                    factor jump*numPixels/dif.     ALERT: if
  243.                    goodIterate < its at any point, we stop
  244.                    cycling */
  245. int     its;
  246. double *pY;
  247. {
  248.     double  MakeAJump();
  249.  
  250.     PreviousAdd = add;
  251.  
  252.     for(;;) {        /* change "add"; if it decreases, don't
  253.                    decrease it too much in one step; but keep
  254.                    re-evaluating while it decreases; decreasing
  255.                    too much causes a problem when two points
  256.                    differ primarily because they are being
  257.                    calculated mod something */
  258.         dif = MakeAJump(its, frac + add, pY);
  259.  
  260.         add =.2 * add +.8 * add * numPixels * jump / dif;
  261.  
  262.  
  263.     /* if add is decreased by more than 20%, go through another time to
  264.        refine add */
  265.         if(add <= 0.5 * PreviousAdd) {
  266.             add = 0.5 * PreviousAdd;
  267.             PreviousAdd = add;
  268.         }
  269.         else
  270.             break;
  271. /*        checkForInterrupt(); seems possibly unnecessary */
  272.     }
  273.     if(add > 2.* PreviousAdd) {/* this multiplier should be at least the
  274.                    reciprocal of multiplier of PreviousAdd */
  275.         add = 2.* PreviousAdd;/* don't increase add too fast */
  276.     }
  277. }
  278.  
  279.  
  280. double  MakeAJump(times, fraction, pY)
  281.                 /* The process is iterated "times" times and
  282.                    compares how far the new y[] is in the plane
  283.                    from the vector pY[], and returns that
  284.                    distance(unless the distance is too small;
  285.                    then it returns .00...01 ); */
  286. int     times;
  287. double  fraction,
  288.        *pY;
  289. {
  290.     double  difX,
  291.             difY,
  292.             ret;
  293.  
  294.     setY(fraction);
  295.  
  296.     process(times);
  297.  
  298.     difX = y[X_coord] - pY[X_coord];
  299.     difY = y[Y_coord] - pY[Y_coord];
  300.     ret = sqrt(difX * difX + difY * difY) / diagonal;
  301.     if(ret <.000001)
  302.         ret =.000001;
  303.     return(ret);
  304. }
  305.  
  306. process(times)            /* the inner loop (which itself iterates
  307.                    multipier times) is iterated up to "iter"
  308.                    times and in fact = "iter" times unless the
  309.                    trajectory leaves the box, in which case
  310.                    iterating is terminated. The position y[] at
  311.                    exit time is its  value after "times"
  312.                    iterates */
  313. int     times;
  314. {
  315.     int     i,
  316.             j;
  317.  
  318.     for(i = 1; i <= times; i++) {
  319.         for(j = 0; j < multiplier; j++)
  320.             (*iteratee) ();
  321.     }
  322. }
  323.  
  324.  
  325.  
  326. int     isInBox(xval, yval)    /* This routine uses the same arithmetic
  327.                    procedures as plot to see if
  328.                    (y[X_coord],y[Y_coord]) is within the window
  329.                    of ScreenSection and if "diameter" is >= 0
  330.                    and it is not within the window, the routine
  331.                    checks to see if it is within
  332.                    "diameters[ScrnSec]" times the width and
  333.                    within "diameters[ScrnSec] times the height
  334.                    of the screen. If it is not close enough to
  335.                    the screen based on this criterion, the
  336.                    program returns NO, otherwise YES */
  337. double  xval,
  338.         yval;
  339. {
  340.  
  341.     if((X_upper - xval) * (X_lower - xval) > ddx
  342.             || (Y_upper - yval) * (Y_lower - yval) > ddy) {
  343.         return(NO);
  344.     }
  345.     return(YES);
  346. }
  347. int     isInCRT(xval, yval)    /* This routine uses the same arithmetic
  348.                    procedures as plot to see if
  349.                    (y[X_coord],y[Y_coord]) is within the window
  350.                    of ScreenSection. If it is not in the
  351.                    screen, the program returns NO, otherwise
  352.                    YES */
  353. double  xval,
  354.         yval;
  355. {
  356.  
  357.     if((X_upper - xval) * (X_lower - xval) > 0
  358.             || (Y_upper - yval) * (Y_lower - yval) > 0) {
  359.         return(NO);
  360.     }
  361.     return(YES);
  362. }
  363.  
  364. boxConstants() {        /* this defines constants needed by isInBox()
  365.                    for checking if((X_upper-xval)*(X_lower-xval) >
  366.                    ddx || (Y_upper-yval)*(Y_lower-yval) > ddy  )
  367.                    the constants ddx and ddy must be reset if
  368.                    the scale or box is changed */
  369.     double  d,
  370.             dd,
  371.             xdif,
  372.             ydif;
  373.  
  374.     xdif = X_upper - X_lower;
  375.     ydif = Y_upper - Y_lower;
  376.     oneOverX = 1 / (X_upper - X_lower);
  377.     oneOverY = 1 / (Y_upper - Y_lower);
  378.     d = diameters[ScrnSec];
  379.     dd = d * (1 + d);
  380.     ddx = dd * xdif * xdif;
  381.     ddy = dd * ydif * ydif;
  382. }
  383.  
  384.  
  385.  
  386. checkForInterrupt() {
  387. #ifndef MAINFRAME
  388.     if(ChkKeyStatus () != 0 || cycle > 0) {
  389.         Interrupt();
  390.         boxConstants();/* this defines constants needed by isInBox()
  391.                    for checking if y[] is in the box, and
  392.                    Interrupt can change the required constants 
  393.                 */
  394.     }
  395. #endif
  396. }
  397.  
  398.  
  399. plotSeg() {            /* plots one or more points depending on
  400.                    numPixels and modFlag */
  401.     double  Xp0,
  402.             Yp0;
  403.     double  jumpSize(), js;
  404.     int     isInCRT(), nPix = numPixels;
  405.     int    i;
  406.  
  407.     if(goodIterate == iter - 1 || oldGoodIt == iter - 1)
  408.         nPix = 10 * nPix;
  409.  
  410.     js = jumpSize();
  411.     if(js >.05 + numPixels *.001)
  412.         return;
  413.  
  414.     if(cross0flag == 1) {    /* this is for preparing to plot a cross */
  415.         plotX = y[X_coord];
  416.         plotY = y[Y_coord];
  417.     }
  418.     if(nPix > 1) {        /* this segment can screw up if either or both
  419.                    points were altered by moduloAB() */
  420.         Xp0 = y[X_coord];
  421.         Yp0 = y[Y_coord];
  422.     /* modFlag == YES means that a number has been changed by the
  423.        moduloAB() */
  424.         for(i = 1; i <= nPix; i++) {
  425.             modFlag = NO;
  426.             plotX = (Xp0 * i + oldX * (nPix - i)) / nPix;
  427.             plotY = (Yp0 * i + oldY * (nPix - i)) / nPix;
  428.  
  429.             if(isInCRT (plotX, plotY) == YES)
  430.                 plot(plotX, plotY);
  431.         }
  432.     }
  433.     else {            /* if numPixels is not greater than 1 */
  434.         modFlag = NO;
  435.         if(isInCRT (y[X_coord], y[Y_coord]) == YES)
  436.             plot(y[X_coord], y[Y_coord]);
  437.     }
  438. }
  439.  
  440.  
  441.  
  442. getGoodIt() {            /* sets y according to frac and iterates as
  443.                    many times up to iter as possible while
  444.                    staying in the box; at the end, y[] is this
  445.                    final image point and goodIterate is the
  446.                    number of iterates */
  447.     double  goodY[EQUATIONS];
  448.     int     isInBox();
  449.     int     i,
  450.             j;
  451.  
  452.     goodIterate = 0;    /* this was = 1 in a buggy earlier version
  453.                    9/7/88 */
  454.  
  455.     setY(frac);
  456.     store(goodY, y, dim);
  457.     for(i = 1; i <= iter; i++) {
  458.         for(j = 0; j < multiplier; j++)
  459.             (*iteratee) ();
  460.  
  461.         if(isInBox (y[X_coord], y[Y_coord]) == YES) {
  462.             store(goodY, y, dim);
  463.             goodIterate = i;
  464.         }
  465.         else {
  466.             store(y, goodY, dim);
  467.             return;
  468.         }
  469.     }
  470. }
  471.  
  472. setY(fraction)
  473. double  fraction;
  474. {
  475.     int     i;
  476.  
  477.     for(i = zeroth; i < zeroth + vec_dim; i++)
  478.         y[i] = ya[i] + fraction * (yb[i] - ya[i]);
  479. /***    y[X_coord] = line_x1 + fraction*(line_x2 - line_x1);
  480.     y[Y_coord] = line_y1 + fraction*(line_y2 - line_y1);   ***/
  481. }
  482.  
  483. double  getAdd(add, pY)    /* returns estimate for add */
  484. double  add,
  485.        *pY;
  486. {
  487.     int     n;
  488.     double  MakeAJump();
  489.  
  490.     getGoodIt();
  491.  /* now estimate add and then refine the estimate */
  492.     for(n = 0; n < 2; n++) {
  493.         dif = MakeAJump(goodIterate, frac + add, pY);
  494.                 /* should be about equal to jump*numPixels */
  495.         add = add * numPixels * USTEP * diagonal / dif;
  496.     }
  497.     return(add);
  498. }
  499. #ifdef OBSOLETE
  500. set_line(xshift1, xshift2, yshift1, yshift2)
  501.                 /* used for drawing lines and unstable
  502.                    manifolds */
  503. double  xshift1,
  504.         xshift2,
  505.         yshift1,
  506.         yshift2;
  507. {
  508.  
  509.     line_x1 = y[eqn0 + X_coord] + xshift1;
  510.     line_x2 = y[eqn0 + X_coord] + xshift2;
  511.     line_y1 = y[eqn0 + Y_coord] + yshift1;
  512.     line_y2 = y[eqn0 + Y_coord] + yshift2;
  513. }
  514. #endif
  515.  
  516. LineForU(iterMultiplier)    /* this routine sets the ends of a line from ya
  517.                    to yb for computing an unstable manifold
  518.                    assuming that period NP in Newton's Method
  519.                    divides the period of the orbit and y1[] is
  520.                    at a point of the orbit */
  521. int     iterMultiplier;
  522. {
  523.     double  diagonal,
  524.             distance(), ratio;
  525.     double  yy[EQUATIONS],
  526.             difX,
  527.             difY,
  528.             dif;
  529.     int     i,
  530.             J;
  531.  
  532.     diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
  533.             + (Y_upper - Y_lower) * (Y_upper - Y_lower));
  534. /*******
  535.     if(printer >= 3)
  536.             PRINT 
  537. "initial:y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
  538. ,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
  539. *********/
  540. /* the following is to set the time coordinates of ya and yb so they will be 
  541.         equal; other coordinates copied are now irrelevant */
  542.     store(ya, y, zeroth);
  543.     store(yb, y, zeroth);
  544.  
  545. /* iterate until y is more than dif/(*USTEP*diagonal) from y1 */
  546.     J = 0;            /* counts iterates to prevent loop from never
  547.                    ending */
  548.     do {
  549.         store(yy, y, dim);
  550.  
  551.         for(i = 0; i < iterMultiplier; i++)
  552.             (*iteratee) ();
  553.  
  554.         difX = y[X_coord] - y[eqn0 + X_coord];
  555.         difY = y[Y_coord] - y[eqn0 + Y_coord];
  556.         dif = sqrt(difX * difX + difY * difY);
  557.         ratio = dif / (numPixels * USTEP * diagonal);
  558.     /* desired ratio = 1 */
  559.         if(printer >= 1 && J == 0 && ratio > 5.)
  560.                 /* ratio should start out small */
  561.             warning1 ();
  562. /*********        if(printer >= 3)
  563.         {
  564.             PRINT 
  565.     "J = %d:  y1(x) = %11.11lf, y() = %11.11lf, dif=%11.11lf\n"
  566.     ,J,y[eqn0+X_coord],y[X_coord],
  567.     y[X_coord]-y[eqn0+X_coord]);
  568.             PRINT 
  569.     "iteration %d: length = %lf length/desired len.= %lf difX=%lf \n"
  570.     ,J,dif,ratio,difX);
  571.         }                    *********/
  572.     } while(++J < 100 && ratio < 1.);
  573.     if(J == 100)
  574.         PRINT
  575.             "Something is wrong. The point may not be unstable \n");
  576.     if(printer >= 2)
  577.         PRINT "diagonal = %lf     ratio = %lf \n", diagonal, ratio);
  578.  
  579. /* now rescale yy so that is iterates to a point dif/(USTEP*diagonal) from y1*/
  580.     if(printer >= 3)
  581.         PRINT "before first rescaling, ratio is now %lf\n", ratio);
  582.     if(ratio > 0) {
  583.         for(i = zeroth; i < lyapzero; i++)
  584.             yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
  585.         store(y, yy, dim);
  586.         for(i = 0; i < iterMultiplier; i++)
  587.             (*iteratee) ();
  588.         difX = y[X_coord] - y[eqn0 + X_coord];
  589.         difY = y[Y_coord] - y[eqn0 + Y_coord];
  590.     }
  591.     dif = sqrt(difX * difX + difY * difY);
  592.     ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
  593.     if(printer >= 3)
  594.         PRINT "after first rescaling, ratio is now %lf\n", ratio);
  595.  
  596. /* now rescale again in an effort to get closer to the desired length */
  597.     for(i = zeroth; i < lyapzero; i++)
  598.         yy[i] = y[eqn0 + i] + (yy[i] - y[eqn0 + i]) / ratio;
  599.     store(y, yy, dim);
  600.     for(i = 0; i < iterMultiplier; i++)
  601.         (*iteratee) ();
  602.     difX = y[X_coord] - y[eqn0 + X_coord];
  603.     difY = y[Y_coord] - y[eqn0 + Y_coord];
  604.     dif = sqrt(difX * difX + difY * difY);
  605.     dif = sqrt(difX * difX + difY * difY);
  606.     ratio = dif / (numPixels * USTEP * diagonal);/* desired ratio = 1 */
  607.     if(printer >= 3)
  608.         PRINT "final ratio is now %lf; about 1 is target \n", ratio);
  609. /*******
  610.     if(printer >= 3)
  611.             PRINT 
  612.     "  J = %d:y1(x) = %11.11lf, y(x) = %11.11lf, dif=%11.11lf\n"
  613.         ,J,y[eqn0+X_coord],y[X_coord],y[X_coord]-y[eqn0+X_coord]);
  614.     if(printer >= 2)
  615.         PRINT "diagonal = %lf final ratio = %lf \n",diagonal, ratio);
  616. ********/
  617. /* now set ends of segment to be ya and yb */
  618.     for(i = zeroth; i < zeroth + vec_dim; i++)
  619.         ya[i] = yy[i];
  620.     for(i = zeroth; i < zeroth + vec_dim; i++)
  621.         yb[i] = y[i];
  622.  
  623.     line_x1 = yy[X_coord];
  624.     line_x2 = y[X_coord];
  625.     line_y1 = yy[Y_coord];
  626.     line_y2 = y[Y_coord];
  627.     connect2 (y[eqn0 + X_coord], y[eqn0 + Y_coord], line_x1, line_y1);
  628. /*    connect2(line_x1,line_y1,line_x2,line_y2);*/
  629.     if(printer >= 2) {
  630.         PRINT
  631.             "Now we compute images of ya-yb = (%lf,%lf)to(%lf,%lf) \n"
  632.             ,line_x1, line_y1, line_x2, line_y2);
  633.         pause(3.);
  634.     }
  635. }
  636.  
  637. LineForB(y1, y2)
  638. double *y1,
  639.        *y2;
  640. {
  641.     diagonal = sqrt((X_upper - X_lower) * (X_upper - X_lower)
  642.             + (Y_upper - Y_lower) * (Y_upper - Y_lower));
  643.     line_x1 = y1[X_coord];
  644.     line_x2 = y2[X_coord];
  645.     line_y1 = y1[Y_coord];
  646.     line_y2 = y2[Y_coord];
  647.     stoplines();         /* for stopping drawing connected lines */
  648. #ifdef X11
  649.     connectp(line_x1, line_y1);
  650. #else
  651.     x_c = line_x1;
  652.     y_c = line_y1;
  653.     plot(x_c, y_c);
  654. #endif
  655.     connectp(line_x2, line_y2);
  656.     return;
  657. }
  658. warning1 () {            /* for lineForU */
  659.     {
  660.         PRINT
  661.             "/n********************************** WARNING ******************************\n"
  662.             );
  663.         PRINT
  664.             " THE FIRST STEP IS BIG. PERHAPS y1 IS NOT AT A PERIODIC ORBIT OF PERIOD NP\n"
  665.             );
  666.         PRINT
  667.             "    PROGRAM WILL ATTEMPT TO CARRY OUT COMMAND\n");
  668.  
  669.         PRINT
  670.             "********************************** WARNING ******************************\n\n"
  671.             );
  672.     }
  673. }
  674.